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Abstract 

Multisite protein modification is a ubiquitous mechanism utilized by cells to control protein 
functions. We have recently proposed a dynamical description of multisite protein modification 
which embodies all the essential features of the process (E. Milotti, A. Del Fabbro, C. Dalla 
Pellegrina, and R. Chignola, Physica A, in press), and we have used this model to analyze the 
stability and the time-scales of this mechanism. The same model can be used to understand 
how the system responds to stimuli: here we show that it displays frequency-dependent dynamical 
hysteresis. This behavior closely parallels - with the due differences - what is observed in magnetic 
systems. By selecting model parameters that span the known biological ranges, we find that the 
frequency-dependent features cover the band of the observed oscillations of molecular intracellular 
signals, and this suggests that this mechanism may have an important role in cellular information 
processing. 
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Multisite protein phosphorylation, nitrosylation, methylation, etc., are very common 
mechanisms in cells, which use these processes to control protein functions and thus 
the flow of information in biochemical networks. We have recently introduced a dynamical 
model of multisite protein modiflcation which describes many properties of these ubiquitous 
mechanisms ||2;]. The basic idea is that a molecule B can dock to anyone of several sites 
on protein A: this produces structural changes in the protein, and eventually, when enough 
sites are occupied by B, the protein undergoes a conformational switch which modifles its 
chemical activity or causes the release of some other substance. In this way it is possible 
to model threshold processes in cells, and the thresholds turn out to be stable and tunable 
over many orders of magnitude. This mechanism is actually quite complex and in [2] we use 
simplifying assumptions that lead to the following differential system for the concentrations: 
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where the states An are the forms of A with exactly n modifled sites out of a total of 
sites, and k+ and k- are the on-off rates. This is a non-dissipative differential system, and 
in [2] we have studied its dynamical properties, we have extended it to include the modifled 
form of A produced by threshold crossing, and flnally found that the bimolecular attach- 
ment/detachment process produces a dynamics that is the same as the classical allosteric 
effect. In this paper we return to the unmodifled system ([T]), and turn our attention to other 
features of the process. 

In (2! we assumed that the concentration of B behaves quasi-statically, so that the fol- 
lowing conservation equation also holds: 

N 

Y,n[An] + [B] = [B]o (2) 

n=l 

where [-B]o is the total concentration of B. Equation ([2]) means that B can either be found 
free in the solution or bound to A, and that each An counts n times as much, because n 



molecules B are bound to it. However equation (j2]) is not adequate to describe a situation 
in which B is driven by an external process. At first sight it may seem that the differential 
equation 

d[B]_ ^ [A^] d[B]o 

n=l 

where d[B]Q/dt is a driving term, solves the problem. However this is not so, because this 
equation can easily push [B] towards unphysical negative values. The problem is that B 
is partly bound to A, and one might drive the total amount of B down to small values so 
fast that the bound B does not have the time to unbind. Biology provides an easy way out: 
the concentration of many proteins is commonly modulated by the interplay of production 
and destruction; production is triggered by genetic expression, while destruction may be 
activated by the expression of a molecule that finally leads to the ubiquitination of the 
target protein, and its degradation by proteasomes. It is easy to see that if the ubiquitination 
process followed by proteasome action is fast enough, the details of the destruction process 
can be neglected, and the change of [B] due to simultaneous production (P) and destruction 
(D) can be described by the equation 

d[B] 
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where VB(t) is the production rate and ku is the decay (destruction) constant associated to 
ubiquitination. When we revise ([3]) to include (jll), we find 

n=l 

Here we use the production rate vsit) as the driving term that modulates the concentration 
of B: we let vsit) change periodically 

VB{t) = VBfi[l + m sin {2iTh't)] (6) 

where m is the modulation index and u is the modulation frequency of the production 
process. This modulation only roughly approximates the actual biological ups and downs 
in enzyme production, but has the obvious advantage that the differential system is excited 
by a single Fourier component, and thus the nonlinearities that are observed in the response 
are unambiguous properties of the differential system. The differential system cannot be 
solved analytically, and we have to resort to numerical integration methods, and explore 



the parameter space as well as we can. On the whole there are eight parameters, i.e., the 
total number of modification sites A^, the total amount of protein A, the on-off rates fc+ and 
k_, the production rate vb,o along with the modulation index m and frequency u, and the 
decay rate kjj. At first sight the exploration of such a large parameter space might seem 
to be a daunting task, but fortunately we can limit ourselves to the restricted ranges of 
the actual biological systems, and we use as starting points some of the values associated 
with the multisite phosphorylation of the retinoblastoma protein (Rb) (as in 
N = 16, [A]o = J2n=i[^n] = lOyuM, k_ = 1 Hz, k+ = 10^ Hz M^^. Using reasonable 
values for ubiquitin concentration and for the forward rate in the Michaelis-Menten reaction 
for ubiquitination, we find that kjj is in the range 1-10000 Hz. From equation we see 
that under stationary conditions, vb = ku[B], and since we know 3 that B has a critical 
breakpoint at [B] = N[A]o, we take vb in the range {0,2Nku[A]Q). The modulation index 
m can only range from to 1, and here we take the full amplitude m = 1. Finally, from 
our previous knowledge of the eigenvalue distribution for the linearized system studied in 
B, we know that the relaxation times in the multisite phosphorylation of the Rb protein 
are approximately in the range 1 ms - 100 s, and therefore - after we have fixed all the 
other parameters - we sweep this range with uniform logarithmic spacing. In the simplified 
situation that we have chosen to highlight the dynamical properties of multisite modification, 
it is natural to monitor the average number of modified sites {n) = Yln=i^l^n]/[A]Q; using 
equation we see that the the static concentration of B should be 

[B]eg = '^[l + msm{27Tiyt)] (7) 



and correspondingly, we expect to find the average number of modified sites 
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Figure [1] shows that this expectation is fulfilled only at very low frequencies: for frequencies 
higher than a few mHz, the curve ([8]) opens up into a loop that is reminiscent of the hysteresis 
loop in magnetic systems. Here the loop cannot be symmetric about the origin, because 
both concentration and the average number of modified sites are bound to be non-negative. 



however, just like in magnetic systems 



J] the loop shape changes and at the same time 



the centroid of the loop migrates to a different position at high frequency. We can define a 
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dynamic order parameter Q as in [J]: 



a 



Q = ^ ^ {n)dt (9) 



which is just (n) averaged over the hysteresis loop, however here it is also useful to consider 
the vertical range R = {n)max — {n)min- Figures [2] and figure [3] show the loop area, the order 
parameter Q, and the vertical range R vs. frequency, for the standard choice of parameters 
listed above. Figure [2] shows that the area behaves roughly as a log-normal function, and 
therefore the high-frequency tail has an approximate power-law dependence jsj. This kind 
of behavior does not change appreciably when the other adjustable parameters are varied 



as explained above. The loop area in magnetic systems in known to follow a scaling law {4]: 
here we find that the loop area obeys the approximate scaling law 

^M«Cexp {-!!«!} (10) 

where the coefficients C, vq, ctq depend on the specific simulation parameters, and therefore 
the tail of the distribution has an approximate 1/v dependence. The coefficients also seem 
to follow simple scaling laws that depend on vb^'- 
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CKo) = Coexp|-^^^^^^^^| (11) 
ooivBfi) = a + bv^^o (12) 
I/O = d + evB^o (13) 

The new coefficients Co, vc, crc, cl, b, w, d, e and r depend on the remaining system 
parameters N, [A]o, k-, and ku. Figure H] shows the superposition of several plots 
of rescaled loop area obtained in different numerical integration runs: the superposition is 
reasonably good, although more work shall be needed to obtain a better scaling law and to 
fix the association between C, uq, cto and the simulation parameters. 

The vertical range R in figure [3] shows that the system behaves as a low pass filter with an 
extended transition region, and since this behavior emerges in a biological setting it is natural 
to wonder whether it may have a deeper meaning. We conjecture that it could be used in 
some biochemical circuits as a kind of slope filter, i.e., a low-pass filter used in frequency 
modulation decoding that converts frequency modulation to amplitude modulation (see, e.g., 
. Figure [5] shows the result of a calculation with a frequency-modulated B production 
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rate 

VB(t) = vb,o {1 + sin [Inuct + msin {2nh'Mt)]} (14) 

where is a fixed frequency, um is tlie modulation frequency, and m is again a modu- 
lation index. Although the system parameters have not been optimized for this purpose, 
we see that indeed a frequency-modulated, fixed-amplitude input produces an amplitude- 
modulated output. 

The importa^nce of hysteresis in enzyme kinetics was first stressed by C. Frieden almost 
40 years ago [7!|: here we have shown how dynamical hysteresis and related effects arise 
naturally in a ubiquitous biological mechanism, multisite protein modification. Hysteresis 
shows up in many biomolecular pathways, and it has been observed experimentally in reac- 
;ions that involve proteins with multiple phosphorylation sites, e.g., the phosphatase Cdc25 
8j. Phosphatase Cdc25 regulates cell division, and it is believed that hysteresis underlies 
the irreversibility of the cell-cycle transition into and out of mitosis. We also note that 
information stored in environmental molecules can propagate in the cell through the oscil- 
lations of various types of intracellular molecules (such as c-AMP, ATP, NADPH, Ca^^) 
and that various important cell functions (e.g., insulin secretion, immune cell activation, 
neuron transmission, gene expression, etc.) are activated by the frequency decoding of this 
information. The Ca^+ oscillations are a paradigmatic instance, and several intracellular 
proteins/enzymes can bind Ca?^ at multiple sites and act as frequency decoders 9j. One 
of these is CaM Kinase II whose enzymatic activity has been experimentally shown to be 
sensitive to the frequency of Ca^^ oscillations Most importantly, the observed periods 
of intracellular Ca^"*" oscillations range between a few seconds to minutes, and match quite 
well the frequency range that we obtain with realistic system parameters (the transition in 
figure [3] spans the range 0.01 Hz - 10 Hz). 
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FIG. 1: Plot of the average number of modified sites (n) vs. [B](,q with system parameters A'^ = 16, 
[A]q = 10/iM, A;_ = 1 Hz, = 10^ Hz M~i. The dotted curve is the static solution The loops 
correspond to different modulation frequencies v: a., v = \ mHz; b. u = mHz; c. v = 100 mHz; 
d. u = 1 Hz. The system follows each loop counterclockwise. 
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FIG. 2: This figure shows the loop area vs. modulation frequency v for = 16, [A\q = 10/iM, 
A;_ = 1 Hz, kj^ = 10^ Hz M~^. The dots represent the individual numerical integrations of the 
differential system. The plot with logarithmic horizontal scale in part b. shows that the curve is 
close to a log-normal shape j^. 
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FIG. 3: Plot a. of the order parameter Q, and b. of the vertical range R vs. modulation frequency 
uiov N = 16, [^]o = 10/iM, k- = l Hz, k+ = 10^ Hz M'^ 
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FIG. 4: Rescaled loop area (arb. units) vs. rescaled frequency x for several different parameter 
sets. The rescaled frequency is x = (v/vq)", where the index a is proportional to I/cq- 
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FIG. 5: Example of amplitude-modulated output obtained from a frequency-modulated fixed- 
amplitude input. In this case the system parameters are as before A'' = 16, [^]o = 10/xM, /j_ = 1 
Hz, k+ = 10^ Hz M-\ and uc = I Hz, um = 0.1 Hz, m = 0.8. 
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